********************************************************************************
** MAIN FIGURES
********************************************************************************

ssc install heatplot coefplot
ssc install palettes, replace
ssc install colrspace


** FIGURE 1 : Estimates across specifications ****
{

	cap mkdir "$figures/Figure1"
	
	foreach model in "lin_ei" "lin_ei_lag" {
		use  "$result_regs/regression_results_wide.dta", clear

		keep if model=="`model'"
		
		foreach yvar in "S" "M" "X" {
			gen ci`yvar'_l = beta`yvar' - 1.96*beta_sd`yvar'
			gen ci`yvar'_h = beta`yvar' + 1.96*beta_sd`yvar'
		}
		gsort -deflator eivar wt bacon control othercontrol 
		
		foreach eivar in "ei" "ei_tot" {
			foreach yvar in "S" {
				preserve
				keep if eivar=="`eivar'"
				drop if (wt=="07" | wt=="10")
				mkmat beta`yvar' ci`yvar'_l ci`yvar'_h, matrix(R)
				mat rownames R = "wt-prod07" ///
					"wt-prod07 lnwage" ///
					"wt-prod07 lnwage exchange" ///
					"wt-prod07 bacon" ///
					"wt-prod07 bacon lnwage" ///
					"wt-prod07 bacon lnwage exchange" ///
					"wt-prod10  " ///
					"wt-prod10  lnwage" ///
					"wt-prod10  lnwage ex" ///
					"wt-prod10 bacon " ///
					"wt-prod10 bacon lnwage" ///
					"wt-prod10 bacon lnwage exchange" ///
					"wt-value07  " ///
					"wt-value07  lnwage" ///
					"wt-value07  lnwage ex" ///
					"wt-value07 bacon " ///
					"wt-value07 bacon lnwage" ///
					"wt-value07 bacon lnwage exchange" ///
					"wt-value10  " ///
					"wt-value10  lnwage" ///
					"wt-value10  lnwage exchange" ///
					"wt-value10 bacon " ///
					"wt-value10 bacon lnwage" ///
					"wt-value10 bacon lnwage exchange" 			
				coefplot (matrix(R[,1]), ci((2 3))), xline(0) graphregion(color(white)) ///
					coeflabels(, notick)
				graph export "$figures/Figure1/coef_`yvar'_`eivar'_`model'.pdf", as(pdf) replace
				restore
			}
			foreach yvar in "M" "X" {
				preserve
				keep if eivar=="`eivar'"
				mkmat beta`yvar' ci`yvar'_l ci`yvar'_h, matrix(R)
				if "`yvar'"=="M" {
					mat rownames R = "wt-imp07  " ///
					"wt-imp07  lnwage" ///
					"wt-imp07  lnwage exchange" ///
					"wt-imp07 bacon " ///
					"wt-imp07 bacon lnwage" ///
					"wt-imp07 bacon lnwage exchange" ///
					"wt-imp10  " ///
					"wt-imp10  lnwage" ///
					"wt-imp10  lnwage exchange" ///
					"wt-imp10 bacon " ///
					"wt-imp10 bacon lnwage" ///
					"wt-imp10 bacon lnwage exchange" ///				
					"wt-prod07  " ///
					"wt-prod07  lnwage" ///
					"wt-prod07  lnwage exchange" ///
					"wt-prod07 bacon " ///
					"wt-prod07 bacon lnwage" ///
					"wt-prod07 bacon lnwage exchange" ///
					"wt-prod10  " ///
					"wt-prod10  lnwage" ///
					"wt-prod10  lnwage exchange" ///
					"wt-prod10 bacon " ///
					"wt-prod10 bacon lnwage" ///
					"wt-prod10 bacon lnwage exchange" ///
					"wt-value07  " ///
					"wt-value07  lnwage" ///
					"wt-value07  lnwage exchange" ///
					"wt-value07 bacon " ///
					"wt-value07 bacon lnwage" ///
					"wt-value07 bacon lnwage exchange" ///
					"wt-value10  " ///
					"wt-value10  lnwage" ///
					"wt-value10  lnwage exchange" ///
					"wt-value10 bacon " ///
					"wt-value10 bacon lnwage" ///
					"wt-value10 bacon lnwage exchange" 	
				}
				if "`yvar'"=="X" {
					mat rownames R = "wt-exp07  " ///
					"wt-exp07  lnwage" ///
					"wt-exp07  lnwage exchange" ///
					"wt-exp07 bacon " ///
					"wt-exp07 bacon lnwage" ///
					"wt-exp07 bacon lnwage exchange" ///
					"wt-exp10  " ///
					"wt-exp10  lnwage" ///
					"wt-exp10  lnwage exchange" ///
					"wt-exp10 bacon " ///
					"wt-exp10 bacon lnwage" ///
					"wt-exp10 bacon lnwage exchange" ///				
					"wt-prod07  " ///
					"wt-prod07  lnwage" ///
					"wt-prod07  lnwage exchange" ///
					"wt-prod07 bacon " ///
					"wt-prod07 bacon lnwage" ///
					"wt-prod07 bacon lnwage exchange" ///
					"wt-prod10  " ///
					"wt-prod10  lnwage" ///
					"wt-prod10  lnwage exchange" ///
					"wt-prod10 bacon " ///
					"wt-prod10 bacon lnwage" ///
					"wt-prod10 bacon lnwage exchange" ///
					"wt-value07  " ///
					"wt-value07  lnwage" ///
					"wt-value07  lnwage exchange" ///
					"wt-value07 bacon " ///
					"wt-value07 bacon lnwage" ///
					"wt-value07 bacon lnwage exchange" ///
					"wt-value10  " ///
					"wt-value10  lnwage" ///
					"wt-value10  lnwage exchange" ///
					"wt-value10 bacon " ///
					"wt-value10 bacon lnwage" ///
					"wt-value10 bacon lnwage exchange" 	
				}
				coefplot (matrix(R[,1]), ci((2 3))), xline(0) graphregion(color(white)) ///
					coeflabels(, notick) 
				graph export "$figures/Figure1/coef_`yvar'_`eivar'_`model'.pdf", as(pdf) replace
				restore
			}	
		}
	}

}


** FIGURE 2 : heatmap for transfer rate and subsidy ****
{
	
	cap mkdir "$figures/Figure2"
	
	* 1. import trasnfer rate results
	use  "$result_regs/regression_results_tr_by_ind.dta", clear	
	merge m:1 spec using "$result_regs/regression_results_wide.dta", ///
		keepusing(eivar bacon wt iv model control othercontrol naics_fe) nogen		
	
	merge m:1 naics using "$result_simul/simulation_dataset.dta", ///
		keepusing(e_d shipments te_2007) keep(3) nogen
	
	* Plot for both direct and comprehensive emissions intensities
	foreach eivar in "ei" "ei_tot" {
	foreach method in "reg" "qreg" {
		preserve
		
		keep if eivar=="`eivar'" | eivar==substr("`eivar'",4,6)
			
		* 4. generate grid for plotting
		local eivar="`eivar'"
		summ e_d, det
		gen ei_grid = (_n-1)*r(max)/1000 if _n <= 1000
		gen te_grid = (_n-1)*1.5/1000 if _n <= 1000
		fillin ei_grid te_grid
				
		* 5. transfer rate plot
		`method' transfer_rate c.e_d##c.te_2007 [iw=shipments]
		gen transfer_rate_grid = ei_grid*_b[e_d] + te_grid*_b[te_2007] ///
												+ te_grid*ei_grid*_b[c.e_d#c.te_2007] + _b[_cons]
		label var transfer_rate_grid "Calibrated Transfer Rate"

		* censoring values for better display of legend
		replace transfer_rate_grid = 1.0 if transfer_rate_grid > 1.0
		replace transfer_rate_grid = 0 if transfer_rate_grid < 0
		
		* 6. subsidy plot
		`method' leakage_rate c.e_d##c.te_2007 [iw=shipments]
		gen subsidy_grid = ei_grid*_b[e_d] + te_grid*_b[te_2007] ///
												+ te_grid*ei_grid*_b[c.e_d#c.te_2007] + _b[_cons]
		label var subsidy_grid "Calibrated Subsidy (tons per $000)"
		
		* censoring values for better display of legend
		replace subsidy_grid = 5.0 if subsidy_grid > 5.0
		replace subsidy_grid = 0 if subsidy_grid < 0
		
		* graph
		replace ei_grid = log(ei_grid)
		gen ln_e_d = log(e_d)
		
		rename transfer_rate_grid TR
		heatplot TR te_grid ei_grid if ei_grid >= 2.5, color(black, intensity(0(.06).9)) cut(-0.05(.1)1.05) ///
				addplot(scatter te_2007 ln_e_d if ln_e_d >= 2.5, msize(tiny) mcolor(gs2)) xlabel(2.5(.5)8.5)  ///
				graphregion(color(white)) aspectratio(1) bins(15) legend(region(lcolor(white))) ///
				xtitle("Logged emissions intensity (EI)") ytitle("Trade exposure (TE)")
		graph export "$figures/Figure2/transfer_rate_`eivar'_`method'.pdf", replace 
		graph export "$figures/Figure2/transfer_rate_`eivar'_`method'.png", replace 
		
		rename subsidy_grid LR
		heatplot LR te_grid ei_grid if ei_grid >= 2.5, color(black, intensity(0(.06).9)) cut(-0.05(.1)1.05) ///
				addplot(scatter te_2007 ln_e_d if ln_e_d >= 2.5, msize(tiny) mcolor(gs2)) xlabel(2.5(.5)8.5)  ///
				graphregion(color(white)) aspectratio(1) bins(15) legend(region(lcolor(white))) ///
				xtitle("Logged emissions intensity (EI)") ytitle("Trade exposure (TE)")
		graph export "$figures/Figure2/subsidy_`eivar'_`method'.pdf", replace
		graph export "$figures/Figure2/subsidy_`eivar'_`method'.png", replace
		
		restore
	}
	}
	
	local i = 1
	foreach f in "subsidy_25" "subsidy_h_p85_25" "subsidy_e_p90_25" {
	
		use  "$result_simul/simulation_allspec_`f'.dta", clear	
		
		* Plot for both direct and comprehensive emissions intensities
		local eivar = "e_d"
		egen median_lr = median(lr_qreg), by(naics eivar)
		keep eivar naics te_2007 subsidy_rate transfer_rate shipments median_lr `eivar'
		
		* generate grid for plotting
		summ `eivar'
		gen ei_grid = (_n-1)*r(max)/1000 if _n <= 1000
		gen te_grid = (_n-1)*1.5/1000 if _n <= 1000
		fillin ei_grid te_grid
		
		gen subsidy = subsidy_rate
		qreg subsidy c.`eivar'##c.te_2007 [iw=shipments]
		gen subsidy_grid = ei_grid*_b[`eivar'] + te_grid*_b[te_2007] ///
												+ te_grid*ei_grid*_b[c.`eivar'#c.te_2007] + _b[_cons]
		label var subsidy_grid "Calibrated TR"
		
		replace subsidy_grid = 1.0 if subsidy_grid > 1.0
		replace subsidy_grid = 0 if subsidy_grid < 0

		if ("`f'"=="subsidy_h_p85_25") {
			gen leakage_risk = transfer_rate * e_d * 25 / 1000.0
			qreg leakage_risk c.`eivar'##c.te_2007 [iw=shipments]
			gen lr_grid = ei_grid*_b[`eivar'] + te_grid*_b[te_2007] ///
												+ te_grid*ei_grid*_b[c.`eivar'#c.te_2007] + _b[_cons]
			_pctile leakage_risk, p(85)
			local limit = r(r1)
			replace subsidy_grid = 0 if lr_grid < `limit'
			replace subsidy_grid = 0.80 if lr_grid >= `limit'
		}
		if ("`f'"=="subsidy_e_p90_25") {
			summ `eivar', det
			local limit = r(p90)
			local limit_low = r(p25)
			summ te_2007, det
			local limit_te = r(p90)
			di `limit_te' 
			di `limit' 
			di `limit_low'
			replace subsidy_grid = 0
			replace subsidy_grid = 0.80 if (ei_grid >=  `limit' ) | (ei_grid >= `limit_low' & te_grid >= `limit_te')
		}

		* graph
		replace ei_grid = log(ei_grid)
		gen ln_e_d = log(e_d)
		if (`i'==1) {
			local nbins = 15
		}
		else {
			local nbins = 60
		}
		summ ei_grid, det
		rename subsidy_grid SR
		heatplot SR te_grid ei_grid if ei_grid >= 2.5, color(black, intensity(0(.06).9)) cut(-0.05(.1)1.05) ///
				addplot(scatter te_2007 ln_e_d if ln_e_d >= 2.5, msize(tiny) mcolor(gs2)) xlabel(2.5(.5)8.5)  ///
				graphregion(color(white)) aspectratio(1) bins(`nbins') legend(region(lcolor(white))) ///
				xtitle("Logged emissions intensity (EI)") ytitle("Trade exposure (TE)")
		graph export "$figures/Figure2/subsidy_sims_`i'.pdf", replace 
		graph export "$figures/Figure2/subsidy_sims_`i'.png", replace 
		rename SR subsidy_grid
		
		local i = `i' + 1
		
	}
	
}
*

** FIGURE 3 Impact of a $10 per Metrix Ton of CO2 Carbon Price ****
{
	
	cap mkdir "$figures/Figure3"
	
	* set tax
	local tax = 25
	
	* start
	use "$result_simul/simulation_allspec_tax_`tax'.dta", clear
	
	* 1) Focus on direct emissions 
	keep if eivar == "ei"
	rename pct_change_p pct_change_p
	
	* 2) Use production elasticities to genereate production response & change in domestic emissions 
	/*collapse (mean) pct_change_p shipments imports exports ///
		(median) beta_prod_p50_=beta_prod beta_imp_p50_=beta_imp beta_exp_p50_=beta_exp beta_con_p50_=beta_con ///
		(p25) beta_prod_p25_=beta_prod beta_imp_p25_=beta_imp beta_exp_p25_=beta_exp beta_con_p25_=beta_con ///
		(p75) beta_prod_p75_=beta_prod beta_imp_p75_=beta_imp beta_exp_p75_=beta_exp beta_con_p75_=beta_con, ///
			by(naics ei)
	
	*/
	
	collapse (mean) pct_change_p shipments imports exports ///
		(median) beta_prod_p50_=beta_prod beta_imp_p50_=beta_imp beta_exp_p50_=beta_exp ///
		(p25) beta_prod_p25_=beta_prod beta_imp_p25_=beta_imp beta_exp_p25_=beta_exp  ///
		(p75) beta_prod_p75_=beta_prod beta_imp_p75_=beta_imp beta_exp_p75_=beta_exp, ///
			by(naics ei)
	
	* scale?
	foreach perc in p25 p50 p75 {
		* production response 
		g pct_change_prod_`perc' = beta_prod_`perc' * pct_change_p 
		g change_prod_`perc' = pct_change_prod_`perc' * shipments // this is DP_i... units = millions of USD?
		
		*g pct_change_con_`perc' = beta_con_`perc' * pct_change_p
		*g change_con_`perc' = pct_change_con_`perc' *  (shipments + imports - exports)
		
	}

	* 3) Generate changes in foreign production
	foreach perc in p25 p50 p75 {
		* import/export changes
		g pct_change_imp_`perc' = beta_imp_`perc' * pct_change_p
		g pct_change_exp_`perc' = beta_exp_`perc' * pct_change_p
		
		g change_imp_`perc' = pct_change_imp_`perc' * imports
		g change_exp_`perc' = pct_change_exp_`perc' * exports
		
	}

	summ pct_change*
	
	gen ei_round = round(ei)
	graph twoway (rspike pct_change_prod_p25 pct_change_prod_p75 ei_round, lcolor(gs6) lwidth(.1)) ///
			(scatter pct_change_prod_p50 ei_round, mcolor(navy*.6) msize(small) msymbol(o) mlwidth(0.01)) ///
			(qfit pct_change_prod_p50 ei_round, color(navy*0.2) lp(dash) lwidth(0.5) ) ///
			, ylabel( , grid gmin gmax glw(thin) glpattern(dot) glc(gs8) tlc(gs8)) //////
			legend(off) ///
			xtitle("Energy Intensity [million Btu/ $ million shipment]") ytitle("Impact of $`tax' per metric ton carbon price") ///
			graphregion(color(white))
	graph export "$figures/Figure3/Figure3_prod.pdf", as(pdf) replace
	
	graph twoway ///
			(rspike pct_change_exp_p25 pct_change_exp_p75 ei_round, lcolor(gs6) lwidth(.1)) ///	
			(scatter pct_change_exp_p50 ei_round, mcolor(dkgreen*.6) msize(small) msymbol(o) mlwidth(0.01)) ///
			(qfit pct_change_exp_p50 ei_round, color(dkgreen*0.2) lp(dash) lwidth(0.5) ) ///
			, ylabel( , grid gmin gmax glw(thin) glpattern(dot) glc(gs8) tlc(gs8)) ///
			legend(off) ///
			xtitle("Energy Intensity [million Btu/ $ million shipment]") ytitle("Impact of $`tax' per metric ton carbon price") ///
			graphregion(color(white))
	graph export "$figures/Figure3/Figure3_exp.pdf", as(pdf) replace
	
	graph twoway ///
			(rspike pct_change_imp_p25 pct_change_imp_p75 ei_round, lcolor(gs6) lwidth(.1)) ///
			(scatter pct_change_imp_p50 ei_round, mcolor(maroon*.8) msize(small) msymbol(o) mlwidth(0.01) ) ///
			(qfit pct_change_imp_p50 ei_round, color(maroon*.2) lp(dash) lwidth(0.5) ) ///
			, ylabel( , grid gmin gmax glw(thin) glpattern(dot) glc(gs8) tlc(gs8)) ///
			legend(off) ///
			xtitle("Energy Intensity [million Btu/ $ million shipment]") ytitle("Impact of $`tax' per metric ton carbon price") ///
			graphregion(color(white))
	graph export "$figures/Figure3/Figure3_imp.pdf", as(pdf) replace
	
}
*


********************************************************************************
** APPENDIX FIGURES
********************************************************************************

cap mkdir $figures/Appendix

	
** FIGURE D.1 ****
/*
Generated by Figure 2 code above.
*/

** FIGURE D.2 ****
/*
Generated by Figure 3 code above.
*/


** FIGURE E.1: Evolution of Energy Prices ****
{

	use "$buildpath/output/regressions_dataset.dta", clear
	keep price_TIV price_imports price_exports shipments imports exports year
	
	* weighted average
	replace price_TIV 		= price_TIV*shipment
	replace price_imports 	= price_imports*imports
	replace price_exports 	= price_exports*exports
	collapse (sum) price_TIV price_imports price_exports shipments imports exports, by(year)
	
	gen price_foreign 		= (price_imports + price_exports)/(imports + exports)
	replace price_imports 	= price_imports/imports
	replace price_exports 	= price_exports/exports
	replace price_TIV 		= price_TIV/shipments
	
	* graph
	graph twoway (line price_TIV year) ///
				 (line price_foreign year), ///
	legend( order(1 "Domestic Price Index" ///
				  2 "Foreign Price Index") col(1)) ///
	ylabel( , grid gmin gmax glw(thin) glpattern(dot) glc(gs8) tlc(gs8)) ///
	ytick(8(1)20 , grid gmin gmax glw(thin) glpattern(dot) glc(gs8) tlc(gs8)) ///
	xlabel(1995(5)2015 , grid gmin gmax glw(thin) glpattern(dot) glc(gs8) tlc(gs8)) ///
	ytitle("Energy Prices ($/MMBtu)" )  xtitle("Year") graphregion(color(white))
	graph export "$figures/Appendix/Figure_EnergyPrices.pdf", as(pdf) replace

}

** FIGURE E.2 **
/*
Not generated, from other sources.
*/


** FIGURE E.3 and E.4: Direct, Indirect and Foreign emissions ****
{
	
	use "$buildpath/input/shipments/NBER_shipments_naics1997_1958-2011.dta", clear
	keep if year==2007
	rename vship shipments
	keep naics shipments
	replace shipments = shipments / 1000 // converts thousands to millions of USD (2007)
	tempfile temp
	save `temp.dta', replace


	* Read in domestic carbon intensities & keep only relevant vars
	use "$EI/emissionintensities_industry_year.dta", clear
	keep year naics CarbonI_TINV
	keep if year==2007
	drop year

	* Merge on domestic energy intensities & keep only relevant vars
	*merge 1:1 naics using "$MAIN/census_energy/generated_energyprice/energy_intensity_industry_year_v3.dta" // 95 matches
	merge 1:1 naics using "$EI/naics_intensity_comp_2007.dta" // 98 matches
	keep if _m==3
	drop _m
	keep naics CarbonI_TINV industrylabel total_intensity intensity

	* Generate US emissions intensities using carbon and energy intensities
	g us_co2_rate_total_combust = CarbonI_TINV * total_intensity / 1000 // tons per million dollars (2007)
	g us_co2_rate_direct_combust = CarbonI_TINV * intensity /1000  // tons per million dollars (2007)
	
	* Merge on foreign emissions intensities
	merge 1:1 naics using "$exiobase/exiobase_emissions.dta" // 98 matches
	keep if _m==3
	drop _m

	* Drop foreign rates for "all" (ie not combustion only)
	drop exp*all* impexp*all* imp*all*

	* Convert foreign units from kg to tons
	foreach v of varlist exp_* imp_* impexp_* {
		replace `v' = `v' / 1000
	}
	
	* Include shipments for weights
	merge 1:1 naics using "`temp.dta'"	
	drop _m
	g baseline_emissions = us_co2_rate_total_combust * shipments

	* US - direct vs. total
	gr twoway (scatter us_co2_rate_total_combust us_co2_rate_direct_combust [w=baseline_emissions], ///
		msymbol(circle_hollow) ///
		xtitle("Direct Emissions Intensity") ytitle("Total Emissions Intensity") ///
		yscale(r(0 5000)) xscale(r(0 5000))) /// 
		(line us_co2_rate_total_combust us_co2_rate_total_combust, legend(off)), ///
		title("Domestic Emissions Intensities: Direct versus Total") note("Units = tons/million USD (2007). Marker size = 2007 baseline emissions.")
	graph export "$figures/Appendix/gr_domestic_compare.png", as(png) replace
	
	* US - direct vs. int'l
	gr twoway (scatter impexp_co2_rate_direct_cmbst_ys us_co2_rate_direct_combust if naics!=325311 [w=baseline_emissions], ///
		msymbol(circle_hollow) ///
		xtitle("Domestic emissions intensity") ytitle("Trade transaction-weighted foreign emissions intensity") ///
		yscale(r(0 5000)) xscale(r(0 5000))) /// 
		(line us_co2_rate_total_combust us_co2_rate_total_combust, legend(off)), ///
		title("Domestic versus Foreign Emissions Intensities") note("Units = tons/million USD (2007). Marker size = 2007 baseline emissions.")
	graph export "$figures/Appendix/gr_foreign_compare.png", as(png) replace		
	
}
*





